* A program that computes the Thomsen Estimator
* See Thomsen (1987) for more detail. 
* Program written by Won-ho Park, University of Michigan
* Copyright Won-ho Park, Aug. 2002 

* The input argument should look like this: 
* thomsen01 var1 var2 var3 var4
* where var1 and var3 are the party votes of interest
* and var2 and var4 are the total number of votes.  
* The variables should all be integers, that is, they should be 
* vote counts, not fractions. 
*
* Also note that the dependent variable (Election 2) is var3
* above.  Thus, if one wants to run a model other than a 
* voter transfer model, than one should specify the independent 
* variable as var3 above. 


*capture discard
*sca drop _all
*mat drop _all
program define ecol4, rclass byable(recall)
  version 7.0
  syntax varlist(min=4 max=4 numeric) [if] [in]
  marksample touse , strok
  tokenize `varlist'
  confirm var `1' `2' `3' `4'
  local vote1 `1' 
  local vote2 `2' 
  local vote3 `3' 
  local vote4 `4' 
**********************************
* Declare Temp variables, macros, matrices
  tempname natsum nat_v1 nat_v2 
  tempname corr n core Loyalty zhi zlow rhi rlow loy_hi loy_lo
  tempname d_core m_natv2 defect Defect de_lo de_hi output props pops

***********************************
* Get the national proportions: 
  quietly  tabstat `vote1' `vote2' `vote3' `vote4' if `touse', s(sum) save
  mat `natsum'=r(StatTot)
  scalar `nat_v1'= `natsum'[1,1]/`natsum'[1,2]
  scalar `nat_v2'= `natsum'[1,3]/`natsum'[1,4]
* nat_v1 and nat_v2 are national proportions. 

**************************************
* Non-linear Tranformation
  tempvar x1 x2 px1 px2 


* x's are party % returns, px's are non-linear transformations like probit or logit
gen `x1'=`vote1'/`vote2' if `touse' 
* Just left it here to check the missing values.
quietly  gen `x2'=`vote3'/`vote4' if `touse'
quietly  gen `px1'=invnorm(`x1')
quietly  gen `px2'=invnorm(`x2')

  quietly corr `px1' `px2' [aw=`vote4'] if `touse'
  scalar `corr' = r(rho)
  scalar `n' = r(N)

* Now do the final integration
sca `core'    = binorm(invnorm(`nat_v1'),invnorm(`nat_v2'), `corr')
sca `Loyalty' = `core'/`nat_v1'
* To derive the standard deviations:
sca `zhi'   = log((1+`corr')/(1-`corr'))/2 + 1.96/sqrt(`n'-3)
sca `zlow'  = log((1+`corr')/(1-`corr'))/2 - 1.96/sqrt(`n'-3)
sca `rhi'   = (exp(2*`zhi' )-1)/(exp(2*`zhi' )+1)
sca `rlow'  = (exp(2*`zlow')-1)/(exp(2*`zlow')+1)
sca `loy_hi'= binorm(invnorm(`nat_v1'),invnorm(`nat_v2'), `rhi' )/`nat_v1'
sca `loy_lo'= binorm(invnorm(`nat_v1'),invnorm(`nat_v2'), `rlow')/`nat_v1'

* The Defection Part
mat `d_core' = `nat_v1'*(`Loyalty',`loy_lo',`loy_hi')
mat `m_natv2'= (`nat_v2', `nat_v2', `nat_v2')
mat `defect' = (`m_natv2' - `d_core')/(1-`nat_v1')
sca `Defect' = `defect'[1,1]
sca `de_lo'  = `defect'[1,3]
sca `de_hi'  = `defect'[1,2]

mat `props'   =(`core', `nat_v2'-`core')
mat `props'   = `props' \ (`nat_v1'-`core', 1-`nat_v2'-`nat_v1'+`core')
*mat list `props'
*mat list `natsum'
mat `pops'   = `props' * `natsum'[1,4]
return matrix Props `props'
return matrix Pops `pops'

mat `output' =(`Loyalty', `loy_lo', `loy_hi')
mat `output' = `output' \ `defect'
* Now get the output
di _newline "{it:Thomsen Estimates, Probit Specification}"
di in text "{hline 19}{c TT}{hline 45}"
di "`vote3'""<{c -}""`vote1'"_col(20)"{c |}" _col(25) "Rates" _col(40)"[95% Conf. Interval]"
di "{hline 19}{c +}{hline 45}"
di "{txt}   Loyalty Rate" _col(20)"{c |}{res}" _col(25) %7.5f `Loyalty' _col(40) "["%7.5f `loy_lo' " ,     " %7.5f `loy_hi' "]"
di "{txt} Defection Rate" _col(20)"{c |}{res}" _col(25) %7.5f `Defect'  _col(40) "["%7.5f `de_lo'  " ,     " %7.5f `de_hi'  "]"
di "{txt}    Correlation" _col(20)"{c |}{res}" _col(25) %7.5f `corr'    _col(40) "["%7.5f `rlow'   " ,     " %7.5f `rhi'    "]"
di in text "{hline 19}{c BT}{hline 45}"
di "{txt}n = {res}" `n'
*===================================================================
* Now finally Store the results in r()
return scalar core = `core'
return scalar rho  = `corr'
* return scalar d_core = d_core' 
return matrix votes `natsum' 
return scalar n = `n'
return matrix out `output'

end
